2106.15163v1 [astro-ph.HE] 29 Jun 2021 


. 
. 


arXiv 


DRAFT VERSION JUNE 30, 2021 
Typeset using IATEX twocolumn style in AASTeX63 


Observation of gravitational waves from two neutron star-black hole coalescences 


THE LIGO SCIENTIFIC COLLABORATION, THE VIRGO COLLABORATION, AND THE KAGRA COLLABORATION 


(Dated: June 30, 2021) 


ABSTRACT 


We report the observation of gravitational waves from two compact binary coalescences in LIGO's 
and Virgo's third observing run with properties consistent with neutron star-black hole (NSBH) 
binaries. The two events are named GW200105.162426 and GW200115_042309, abbreviated as 
GW200105 and GW200115; the first was observed by LIGO Livingston and Virgo, and the second 
by all three LIGO- Virgo detectors. The source of GW200105 has component masses Bore Me and 
1.9*03 Mo, whereas the source of GW200115 has component masses 5.7 3:1 Mo and 1.5703 Mo (all 
measurements quoted at the 9096 credible level). The probability that the secondary's mass is below 
the maximal mass of a neutron star is 8996-9696 and 8796-9896, respectively, for GW200105 and 
GW200115, with the ranges arising from different astrophysical assumptions. The source luminosity 
distances are 2801110 Mpc and 3001150 Mpc, respectively. The magnitude of the primary spin of 
GW200105 is less than 0.23 at the 90% credible level, and its orientation is unconstrained. For 
GW200115, the primary spin has a negative spin projection onto the orbital angular momentum at 
88% probability. We are unable to constrain the spin or tidal deformation of the secondary component 
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for either event. We infer an NSBH merger rate density of 45:33 Gpc ? yr ! when assuming that 
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GW200105 and GW200115 are representative of the NSBH population, or 130766 Gpc 3 yr! under 
the assumption of a broader distribution of component masses. 


1. INTRODUCTION 


In January 2020, the LIGO-Virgo detector net- 
work observed gravitational-wave (GW) signals from 
two compact binary inspirals which are consistent 
with neutron star-black hole (NSBH) binaries. These 
represent the first confident observations to date of 
NSBH binaries via any observational means. The two 
events, carrying the full designations GW200105_162426 
and GW200115_042309, and abbreviated henceforth as 
GW200105 and GW200115, were detected on January 
05, 2020 at 16:24:26 UTC and January 15, 2020 at 
04:23:10 UTC, respectively. The coincident detection 
of GW200115 by the three detectors LIGO Hanford, 
LIGO Livingston, and Virgo gives it high confidence 
of being an astrophysical GW event. During the other 
event, GW200105, the LIGO Hanford detector was not 
operational and, owing to the small signal-to-noise ra- 
tio (S/N) in Virgo, it was effectively a single-detector 
event in LIGO Livingston. While quantification of the 
confidence of single-detector events is subject to signifi- 
cant uncertainty, GW200105 stands clearly apart in the 
LIGO Livingston data from any other candidate with 


NSBH-like parameters during the ~11 months of the 
third observing run (O3). 

'The component masses inferred from these binary in- 
spirals provide information on the nature of their com- 
ponents. The primaries have masses of mı = 8.9 12 Mo 
and mı = 5.7*1* Mo, for GW200105 and GW200115, 
respectively, with uncertainties quoted at the 9096 cred- 
ible level. These primary masses are well above the max- 
imum mass of a neutron star (NS; Rhoades & Ruffini 
1974; Abbott et al. 2018a; Shibata et al. 2019; Cromar- 
tie et al. 2019; Farr & Chatziioannou 2020; Nathanail 
et al. 2021; Fonseca et al. 2021) and within the mass- 
range of black holes (BHs) observed electromagnetically 
(Ozel et al. 2010; Farr et al. 2011; Kreidberg et al. 
2012; Miller & Miller 2014) and via GWs (Abbott et al. 
2016a, 2019a, 2021c). The masses of the secondaries are 
ma —1.9*03 Mo and ma —1.5*0 3 Mo, respectively, for 
GW200105 and GW200115, within the mass range of 
known NSs (Antoniadis et al. 2016; Alsing et al. 2018; 
Abbott et al. 20178, 20208). 

Detections of NSBHs have so far remained elusive in 
both electromagnetic (EM) and GW surveys. In the 
past four decades, surveys have identified 19 binary neu- 
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tron star (BNS) systems in the Galaxy (Farrow et al. 
2019; Agazie et al. 2021); however, the discovery of a 
pulsar in an NSBH binary remains a key objective for 
current and future radio observations (Liu et al. 2014; 
Weltman et al. 2020). 

Similarly, GW observations of LIGO and Virgo 
through the first part of the 03 run (03a) have led to 
the identification of 48 binary black hole (BBH) candi- 
dates (Abbott et al. 2019a, 2021c) and two BNS candi- 
dates (Abbott et 81. 2017a, 2020c). Independent anal- 
yses of the public detector data identified additional 
GW candidates (Venumadhav et al. 2019, 2020; Za- 
ckay et al. 2019a,b; Magee et al. 2019; Nitz et al. 2019, 
2020b; Nitz et al. 2021). The absence of NSBH candi- 
dates in LIGO and Virgo's first and second observing 
runs (O1 and O2, respectively) led to the upper limit 
on the local merger rate density of NSBH systems of 
RxsBu < 610 Gpc * yr^! at the 90% credible level (Ab- 
bott et, al. 2019a). 

During O3a, two events were notable as possible 
NSBH candidates. First, GW190426 152155 (Abbott 
et al. 2021c) was identified as a marginal NSBH can- 
didate with false alarm rate (1.4/yr) so high that it 
could also plausibly be a detector noise artifact. Sec- 
ond, GW190814 (Abbott et al. 2020d) may have been 
an NSBH merger. Although GW190814’s secondary 
mass of mo = 2.59*005 Mo likely exceeds the maxi- 
mum mass supported by slowly spinning NSs (Essick 
& Landry 2020; Fattoyev et al. 2020; Tews et al. 2021; 
Godzieba et al. 2021), such as those found in known 
binaries that will coalesce within a Hubble time, the 
secondary could conceivably be an NS spinning near its 
breakup frequency (Essick & Landry 2020; Tews et al. 
2021; Most et al. 2020; Dexheimer et al. 2021). 

'The existence of NSBH systems has long been con- 
jectured. Observations in the Milky Way reveal high- 
mass X-ray binaries composed of a massive star and a 
compact object (Liu et al. 2008; Orosz et al. 2009; Gou 
et al. 2009; Orosz et al. 2011; Gou et al. 2014). Binary 
evolution models show that X-ray binaries with a BH 
component are possible progenitors of NSBH systems 
(Belczynski et al. 2013; Grudzinska et 81. 2015). 

Major uncertainties regarding massive binary evolu- 
tion, such as mass loss, mass transfer and the impact of 
supernova explosions result in a wide range of merger 
rate predictions: 0.1-800 Gpc^? yr^! (Belezynski et al. 
2002; Sipior & Sigurdsson 2002; Belczynski et al. 2006; 
Dominik et al. 2015; Belczynski et al. 2016; Eldridge 
et al. 2017; Mapelli & Giacobbo 2018; Kruckow et al. 
2018; Neijssel et al. 2019; Drozda et al. 2020; Zevin 
et al. 2020; Broekgaarden et al. 2021). Comparable 
NSBH merger rates are predicted from young star clus- 


ters (Ziosi et al. 2014; Rastello et al. 2020) while NSBH 
merger rates in globular and nuclear clusters are pre- 
dicted to be orders of magnitude lower (Clausen et al. 
2013; Ye et al. 2020; Fragione & Banerjee 2020; Hoang 
et al. 2020; Arca Sedda 2020). Measuring the NSBH 
merger rate and properties such as masses and spins is 
crucial in determining formation channels. 

'This letter presents the status of the detectors dur- 
ing times around GW200105 and GW200115 (82), and 
the results of the different searches leading to the de- 
tections (83). We describe the main properties of the 
two events (84) and discuss the nature of the secondary 
components (85). Finally, we present the astrophysical 
implications, including merger rates of this new class of 
GW source (86) and conclude (87). Data products asso- 
ciated with the events reported here, such as calibrated 
strain time series and parameter estimation posterior 
samples, are available through the Gravitational Wave 
Open Science Center (GWOSC) at gw-openscience.org. 


2. DETECTORS AND DATA 


'The two events reported here were observed in the sec- 
ond part of the third observing run (O3b). At the time 
of GW200105, LIGO Livingston had been in a stable op- 
erational state for over 10 hrs, with a sensitivity, quan- 
tified by the angle-averaged BNS inspiral range (Allen 
et al. 2012), of 137 Mpc. Virgo had been in its nom- 
inal state for ~22 hrs, with a BNS range of ~45 Mpc, 
and LIGO Hanford was not operational. At the time 
of GW200115, all three interferometers had been in a 
stable operational state for over 2 hrs. The BNS ranges 
for LIGO Hanford, LIGO Livingston, and Virgo around 
the detection were ~115 Mpc, ~133 Mpc, and “50 Mpc, 
respectively. 

Figure 1 shows time-frequency representations (Chat- 
terji et al. 2004) of the two events. The LIGO Livingston 
data of GW200105 show a track of excess power with 
increasing frequency. For GW200115, no similar tracks 
are visible, as the S/N in each of the detectors is lower 
than that of GW200105 in LIGO Livingston (see Fig. 3). 
Also, light-scattering noise (Soni et al. 2020; Soni et al. 
2021) is visible in Fig. 1 in LIGO Livingston around 
20 Hz. 

'The LIGO and Virgo GW detectors are calibrated us- 
ing radiation pressure from auxiliary lasers at a known 
frequency and amplitude (Acernese et al. 2018; Sun et al. 
2020). In LIGO, the calibration pipeline (Viets et al. 
2018) linearly subtracts the noise from the calibration 
lines and the harmonics of the power mains from the 
data. 

The calibration systematic error and associated uncer- 
tainty of the data in the [20-1024] Hz frequency region 
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Time-frequency representations of the data containing GW200105 (left column) and GW200115 (right column). 


Times are shown relative to the signals’ merger times, January 05, 2020 at 16:24:26 UTC (left) and January 15, 2020 at 04:23:10 
UTC (right). The amplitude scale of each time-frequency tile is normalized by the respective detector's noise amplitude spectral 
density. The LIGO Livingston data of GW200105 show a track of excess power with increasing frequency. In the other panels, no 
similar tracks are visible, as the S/N in each of the detectors is lower and (for GW200115) the signal is longer. For GW200105, 
the LIGO Livingston data are shown after glitch subtraction. For GW200115, light-scattering noise is visible in LIGO Livingston 


below 25 Hz. 


are for the LIGO detectors no larger than 8.696 in am- 
plitude and 5.9? in phase for GW200105 and 8.096 in 
amplitude and 5.5? in phase for GW200115 (6896 cred- 
ible intervals). For Virgo, we use 5.0% in amplitude 
and 4.1? in phase for both events, except for frequen- 
cies 46-51 Hz, where additional calibration systematic 
error arises from online loops set up to damp the main 
power lines at 50 Hz and mechanical resonances of the 
suspensions close to 48 Hz. This effect was introduced 
during detector improvements carried out between O3a 
and O3b and is not accounted for in the signal recon- 
struction process. However, Virgo data in the frequency 
window 46-51 Hz are excluded from the parameter esti- 
mation analyses in 84. 

'To verify that instrumental noise artifacts do not bias 
the analysis of source properties of the observed events, 
we use data quality validation procedures as in previ- 
ous events (Abbott et al. 2016b; Davis et al. 2021), em- 
ploying sensor arrays at LIGO and Virgo to measure 


environmental disturbances that could couple into the 
interferometers (Nguyen et al. 2021). In Virgo, we find 
no evidence of excess power from terrestrial sources for 
both events. For GW200105, we identify light-scattering 
noise in LIGO Livingston below 25 Hz, 3 s before merger. 
As in past detections (Abbott et al. 2021c), we sub- 
tract this noise with the BayesWave algorithm (Cornish 
et al. 2021) and use the cleaned data for the source pa- 
rameter estimation in 84. The top left panel of Fig. 1 
shows the cleaned data. A low-energy feature around 
20 Hz, not overlapping with the time-frequency track 
of GW200105, remains after the glitch subtraction. For 
GW200115, we also identify light scattering in the LIGO 
Livingston data below 25 Hz. Due to the increased dif- 
ficulty of subtracting the long-duration glitching that 
coincides with the time-frequency track of GW200115, 
glitch-subtracted data were not available at the time of 
analysis. Hence, we exclude LIGO Livingston data be- 
low 25 Hz in the analysis in 84. 
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Figure 2. Sky localizations for GW200105 (top) and for 
GW200115 (bottom), in terms of right ascension and declina- 
tion. The thick, solid contours show the 9096 credible regions 
from the low-latency sky localization algorithm BAYESTAR 
(Singer & Price 2016). The shaded patch is the sky map ob- 
tained from the preferred high-spin analysis of 84, with the 
9096 credible regions bounded by the thin dotted contours. 


3. DETECTIONS 


GW200115 is a coincident event, and in 83.1 we de- 
scribe the established procedures to identify it and 
determine its significance. Subsequently, we describe 
the procedures followed for the single-detector event 
GW200105. 


3.1. GW200115 - Multi-detector event 


GW200115 was initially identified by all four low- 
latency matched-filtering pipelines as a possible com- 
pact binary coalescence (CBC) candidate in LIGO Han- 
ford and LIGO Livingston: GsTLAL (Cannon et al. 
2012; Privitera et al. 2014; Messick et al. 2017; Sachdev 
et al. 2019; Hanna et al. 2020), MBTAONLINE (Aubin 
et al. 2021; Adams et al. 2016), PYCBC Live (Dal Can- 
ton et al. 2020; Nitz et al. 2018, 2019, 2017; Usman 
et al. 2016), and SPIIR (Luan et al. 2012; Hooper et al. 
2012; Chu et al. 2020; Guo et al. 2018). As detailed 
in Abbott et al. (2021c), matched-filtering searches use 
banks (Owen & Sathyaprakash 1999; Harry et al. 2009; 
Privitera et al. 2014; Roy et al. 2017, 2019) of modeled 
gravitational waveforms, with the mass range relevant 
to GW200105 and GW200115 covered by the SEOB- 


Table 1. Network S/N recovered by the search pipelines 
in low-latency and in later offline analysis with improved 
calibration and refined data-quality information. An asterisk 
(*) indicates values where Virgo is not included in the S/N 
calculation. 


Event GstTLAL MBTA PvCBC SPIIR 


GW200105 low-latency 13.9 13.3 13.2 13.2 
offline 13.9 13.4 13.1* — 

GW200115 low-latency 11.4 11.4 11.3 11.0 
offline 11.6 11.2 10.8* — 


NRv4_ROM waveform model (Bohé et al. 2017; Pürrer 
2016). The signal in Virgo was not loud enough to fur- 
ther improve the significance of the coincident candidate 
observed by the LIGO detectors, so GW200115 was re- 
ported as a two-detector event in the GRACEDB event 
database named $200115j with a merger time of January 
15, 2020 at 04:23:10 UTC. 

A public Preliminary Gamma-ray burst Coordinates 
Network (GCN) Notice was sent out six minutes af- 
ter the event (LIGO Scientific Collaboration & Virgo 
Collaboration 2020a). The low-latency BAYESTAR 
(Singer & Price 2016) sky map computed from the orig- 
inal trigger, which was produced from MBTAONLINE, 
indicated a possible discrepancy compared to those from 
the other three pipelines. Therefore, the representative 
trigger was manually switched to the one from GSTLAL 
(LIGO Scientific Collaboration & Virgo Collaboration 
2020b). The associated sky map is shown in the bottom 
panel of Fig. 2 with a solid black line with 90% credible 
area of 900 deg?. The GCN Circular reported a prob- 
ability > 99% for the lighter compact object to have a 
mass below 3Mo, derived from a machine-learning anal- 
ysis (Kapadia et al. 2020; Chatterjee et al. 2020). The 
RAVEN pipeline (Urban 2016), which looks for coinci- 
dent gamma ray bursts (GRBs), did not report an asso- 
ciated GRB trigger from Fermi-GBM or the Neil Gehrels 
Swift Observatory within —1 s and 4-5 s of the GW trig- 
ger. A total of 31 GCN Circulars (GCN archive for 
S200115j 2020) reporting follow-up observations were 
published for this event. To date, no EM counterpart 
has been reported associated with GW200115. 

After the low-latency identification of GW200115, the 
extended periods of strain data around the event are 
further analyzed by the detection pipelines in their of- 
fline configurations using improved calibration and re- 
fined data quality information that is not available in low 
latency (Abbott et al. 2018b, 2021c; Davis et al. 2021). 
In this analysis, we also subtract nonstationary noise 
due to the nonlinear coupling of the 60 Hz power mains 


at the LIGO detectors using coupling functions derived 
from machine-learning techniques (Vajente et al. 2020). 

PYCBC uses a template bank constructed with a 
hybrid geometric-random algorithm (Roy et al. 2019, 
2017), while GSTLAL and MBTA use a template bank 
generated by a stochastic placement method (Mukher- 
jee et al. 2021; Babak 2008; Privitera et al. 2014). GST- 
LAL identifies GW200115 with a network S/N of 11.6 
and FAR of <1/(1x 10° yr) using the data from Novem- 
ber 1, 2019 15:00 UTC to January 22, 2020 18:11 UTC. 
The MBTA and PvCBC offline analyses also identified 
the trigger in data from January 13, 2020 10:27 UTC 
to January 22, 2020 18:11 UTC, yielding consistent net- 
work S/Ns as shown in Table 1, as well as estimated 
FARs of 1/(182 yr) and <1/(5.6 x 104 yr), respectively. 
Based on two detectors, these FARs robustly indicate 
the confidence of the detection. 


3.2. GW200105 — Single-detector event 


GW200105 was identified by GSTLAL running in the 
low-latency configuration as a possible CBC candidate 
in LIGO Livingston and Virgo data with network S/N of 
13.9 and merger time January 05, 2020 at 16:24:26 UTC. 
The candidate identified as 520010586 in GRACEDB was 
considered a single-detector event because the S/N in 
Virgo was below the threshold S/N of 4.0. 

'The low-latency FAR of z 24 per year did not pass 
the threshold for sending a GCN alert. Without the 
information provided by temporal coincidence between 
different detectors and consistency in recovered param- 
eters, it is more difficult to obtain a robust estimate 
of the event’s significance. In fact, an alternative con- 
figuration of GSTLAL, running to test several improve- 
ments made in the offline configuration of early O3 (Ab- 
bott et al. 2021c), found the trigger at higher signifi- 
cance. Therefore, validation procedures and subsequent 
early-offline analysis were initiated. The three other low- 
latency search pipelines, MBTAONLINE, PYCBC LIve, 
and SPIIR, also generated triggers with consistent S/Ns 
near the event time of S200105ae (see Table 1). These 
three pipelines were not configured to assign FARs to 
single-detector triggers and therefore did not generate 
automatic alerts for GW200105. 

On January 6, 2020 at 19:39:09 UTC, the initial GCN 
Circular reported a low-latency BAYESTAR skymap 
(solid black contours) in the top panel of Fig. 2 with 90% 
credible region of 7700 deg? and a 1296 chance of for- 
mation of tidally disrupted matter—relevant for a pos- 
sible EM counterpart—as a result of the coalescence, 
as derived from a machine-learning analysis (Kapadia 
et al. 2020; Chatterjee et al. 2020). Using the param- 
eters obtained from low-latency parameter estimation, 
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Figure 3. Colored shading shows the joint S/N-£? 
noise probability density function for LIGO Hanford (LHO), 
LIGO Livingston (LLO), and Virgo, computed from back- 
ground triggers found by GSTLAL in the region of lighter 
binary systems during the entire O3. The red star indicates 
GW200105. At its position, there is no background present, 
indicating that it stands above all of the recorded background 
triggers. For comparison, the triggers of GW200115 and the 
marginal GW190426.152155 are also shown. 


the chance of formation of tidally disrupted material 
was later revised to be negligible (see LIGO Scientific 
Collaboration & Virgo Collaboration 2020c and 85 be- 
low). 

In response to the discovery notice of GW200105, mul- 
tiple observatories carried out follow-up and shared their 
results publicly via 21 GCN circulars (GCN archive for 
820010586 2020). To date, no EM counterpart has been 
reported associated with GW200105. 

After the event's identification in low latency, the 
strain data are further analyzed by the detection 
pipelines in their offline configuration, following pro- 
cedures for single-detector events discussed in Abbott 
et al. (20208). GSTLAL identifies GW200105 with an 
S/NR of 13.9 and FAR of 1/(2.8 yr) in its offline config- 
uration using data from November 1, 2019 15:00 UTC 
up to January 22, 2020 18:11 UTC. For single-detector 
events, like GW200105, the FAR estimate involves ex- 
trapolation, as explained below. The MBTA and Pv- 
CBC offline analyses also identify the trigger in the 
LIGO Livingston data from January 4, 2020 17:09 UTC 
to January 13, 2020 10:27 UTC yielding the consistent 
S/Ns shown in Table 1. The S/N for GW200105 is larger 
than that for GW200115, even though LIGO Hanford 
was not operational. 

We cannot rely on only the S/N of the CBC trigger to 
estimate its significance, i.e. to quantify how often detec- 
tor noise mimics a possible coalescence signal. Because 
detector noise shows significant deviation in the tails 
from the standard Gaussianity assumptions, its prop- 
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erties have to be estimated empirically (Abbott et al. 
2016c). For multi-detector triggers like GW200115, con- 
sistency of the observed CBC trigger among two or more 
detectors forms an integral part of establishing it as 
a confident detection with a low FAR. On the other 
hand, single-detector triggers can be assigned only mod- 
est FAR values due to limited detector observational 
time (Callister et al. 2017; Nitz et al. 2020a). Nonethe- 
less, alternative methods can be used to estimate con- 
fidence in a single-detector event, as was demonstrated 
for GW190425 (Abbott et al. 2020a). 

The signal GW200105 measured in LIGO Livingston 
is distinct from all noise events in the entire 03 obser- 
vation period. Specifically, Fig. 3 shows with colored 
shading the distribution of background noise triggers 
identified by GSTLAL in the region of binary systems 
with a chirp mass less than 4 Mọ. The colored region 
indicates the density of background noise triggers quan- 
tified through the S/N-£? noise probability density func- 
tion for the three detectors, where the autocorrelation £? 
provides a consistency test similar to x? (Messick et al. 
2017). The red star represents GW200105, which lies in 
a region of this plane without noise triggers, indicating 
that this trigger is unique within the entire O3. 

For comparison, Fig. 3 also shows GW200115 and 
the marginal candidate GW190426_152155 separately 
for LIGO Livingston and LIGO Hanford. In general, 
triggers such as GW200115 and GW190426_152155 in 
Fig. 3 would not be separable from the noise by a single 
detector alone, and it is the coincidence between multi- 
ple detectors that raises their significance. On the other 
hand, the coincidence was not available for GW200105, 
yet it clearly stands out from the background in Fig. 3 
as strong GW signals typically do. 

'The uniqueness of GW200105 within the data shown 
leads to a upper bound on the FAR assignment compa- 
rable to the inverse observing time. A stronger bound, 
like the FAR quoted above, must rely on assumptions 
made on the properties of the noise triggers if one had 
collected the observational strain data for a longer pe- 
riod, a process called extrapolation. As a consequence of 
the assumptions made, the extrapolated FAR estimates 
have large uncertainties. Since single-detector FAR es- 
timates rely on the uniqueness within the dataset, the 
values may change and errors may reduce as more data 
are accumulated. At the time of the analysis, nei- 
ther MBTA nor PvCBC had the capability to as- 
sign a significance to single-detector triggers, although 
GW200105 also stands out as a unique trigger in their 
analyses. Based on these considerations, we consider 
this trigger to be astrophysical and include it in the 
analysis in the remainder of this paper. 


4. SOURCE PROPERTIES 


We infer the physical properties of the two GW 
events using a coherent Bayesian analysis following the 
methodology described in Appendix B of Abbott et al. 
(20198). For GW200105, data from LIGO Livingston 
and Virgo are analyzed, whereas for GW200115, data 
from both LIGO detectors and Virgo are used. 

Owing to the different signal durations, we analyze 
32s of data for the higher-mass event GW200105 and 
64s of data for GW200115. All likelihood evaluations 
use a low-frequency cutoff of fiow = 20 Hz, except for 
LIGO Livingston for GW200115, where fij, = 25 Hz 
avoids excess noise localized at low frequencies, as dis- 
cussed in 82. The power spectral density used in the 
likelihood calculations is the median estimate calculated 
with BayesLine (Cornish & Littenberg 2015). 

The parallel Bilby (PBILBY) inference library, together 
with the DYNESTY nested sampling software (Ashton 
et al. 2019; Romero-Shaw et al. 2020; Smith et al. 2020; 
Speagle 2020) is the primary tool used to sample the 
posterior distribution of the sources! parameters and 
perform hypothesis testing. In addition, we use RIFT 
(Lange et al. 2018) for the most computationally expen- 
sive analyses and LALINFERENCE (Veitch et al. 2015) 
for verification. 

We base our main analyses of GW200105 and 
GW200115 on BBH waveform models that include the 
effects of spin-induced orbital precession and higher- 
order multipole GW moments, but do not include tidal 
effects on the secondary. Specifically, we use two sig- 
nal models: IMRPhenomXPHM (Phenom PHM; Prat- 
ten et al. 2020) from the phenomenological family and 
SEOBNRv4PHM (EOBNR PHM; Ossokine et al. 2020) 
from the effective one-body numerical relativity family. 
'The acronym PHM stands for Precessing Higher-order 
multipole Moments. Henceforth, we will use the short- 
ened names for the waveform models. 

In order to quantify the impact of neglecting tidal ef- 
fects, we also analyze GW200105 and GW200115 us- 
ing two NSBH waveform models that include tidal ef- 
fects and assume that spins are aligned with the or- 
bital angular momentum: IMRPhenomNSBH (Phe- 
nom NSBH; Thompson et al. 2020) and SEOB- 
NRv4. ROM. NRTidalv2 NSBH (EOBNR NSBH; Matas 
et al. 2020). We restrict the NSBH analyses to the region 
of applicability of the NSBH models, i.e. x1 < 0.5, X2 < 
0.05 for Phenom NSBH and x, < 0.9, xa < 0.05 for 
EOBNR NSBH. We also perform aligned-spin BBH 
waveform analyses and find good agreement with the 
analyses using NSBH waveform models (see 84.6 below), 
validating the use of BBH waveform models. Specif- 
ically, we use the aligned-spin BBH models IMRPhe- 
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Table 2. Source properties of GW200105 and GW200115. We report the median values with 90% credible intervals. Parameter 
estimates are obtained using the Combined PHM samples. 


GW200105 GW200115 

Low Spin High Spin Low Spin High Spin 

(x2 < 0.05) (x2 < 0.99) (x2 < 0.05) (xa < 0.99) 
Primary mass mi/Mo y 8.9*12 5.9514 5.718 
Secondary mass m2/Mo iE 1.9755 loger 15707 
Mass ratio q 0.214008 022504 0.244931 020017. 
Total QUARE -M MS tto  10:9513 78:15. TAMA 
Chirp mass M/Mo 3.41007 3.415557 24270 2421 00ኛ 


Detector-frame chirp mass (1+ z)M/Mo 


Primary spin magnitude xi 
Effective inspiral spin parameter Xeff 


Effective precession spin parameter xp 


Luminosity distance Dr /Mpc 


3.61915 a8 
—0.01 


280711 


3.619 6.008  2-580Ż0007 2.579 600 


0.09:0:08.. 0.0820.05 (up - tages 
+008 (.ተ Spd Sater te 
0.07*6:06 000920. 0107." 0.211939 
2807115 310 99 =. 300 +150 

0.067 5:5 (unes Hake 


Source redshift z 


nomXAS (Phenom; Pratten et al. 2020a) and SEOB- 
NRv4 (EOBNR; Bohé et al. 2017), which only contain 
dominant quadrupole moments, and IMRPhenomXHM 
(Phenom HM; García-Quirós et al. 2020) and SEOB- 
NRv4HM (EOBNR HM; Cotesta et al. 2018, 2020), 
which contain higher-order moments. 

The secondary objects are probably NSs based on 
mass estimates, as discussed in detail in 85. As in ear- 
lier GW analyses (Abbott et al. 2017a, 2020a), we pro- 
ceed with two different priors on the secondary's spin 
magnitude: a low-spin prior, X2 X 0.05, which captures 
the maximum spin observed in Galactic BNSs that will 
merge within a Hubble time (Burgay et al. 2003), and a 
high-spin prior, x2 < 0.99, which is agnostic about the 
nature of the compact object. The two priors allow us 
to investigate whether the astrophysically relevant sub- 
case of low NS spin leads to differences in the parameter 
estimation for the binaries. All other priors are set as in 
previous analyses (e.g., Abbott et al. 2021c). Through- 
out, we assume a standard flat ACDM cosmology with 
Hubble constant Ho = 67.9kms-! Mpc ! and matter 
density parameter Qm = 0.3065 (Ade et al. 2016). 

For each spin prior, we run our main analyses with 
higher-order multipole moments and precession for both 
waveform families, EOBNR PHM and Phenom PHM. 
The EOBNR PHM model is used in combination with 
RIFT and the Phenom PHM model with PBILBy. The 
parameter estimation results for the individual precess- 
ing waveform models yield results in very good agree- 
ment; the median values typically differ by 1/10 of the 
width of the 90% credible interval. Nevertheless, in or- 
der to alleviate potential biases due to different sam- 


0.06" $0 


plers or waveform models, we combine an equal num- 
ber of samples of each into one data set for each spin 
prior (Abbott et al. 2016d; Ashton & Khan 2020; Ab- 
bott et al. 2020d) and denote these as Combined PHM. 
The quoted parameter estimates in the following sec- 
tions are the Combined PHM high-spin prior analyses. 
In the figures, we emphasize the high-spin prior results. 
The values of the most important parameters of the bi- 
naries are summarized in Table 2, and we will present 
details in the following sections. 


4.1. Masses 


Figure 4 shows the posterior distribution for the com- 
ponent masses of the two binaries. Defining the mass 
parameters such that the heavier mass is the primary 
object, i.e. m4 > m», our analysis shows that GW200105 
is a binary with a mass ratio of q = m2/m, = UP PRI ተ 
with source component masses m, = 8.91 Mo and 
mə = 1.9*03 Mo. Similarly, GW200115 is a binary with 
a mass ratio of q = 265: with source component 
masses mı = 5.7715 Mo and mə = Lor Mo. 

The primary components of GW200105 and 
GW200115 are identified as BHs from their mass mea- 
surements. For GW200115, we find that the prob- 
ability of the primary falling in the lower mass gap 
(3 Mo ፎ mi ፎ 5 Mo; Bailyn et al. 1998; Orosz et al. 
1998) is 30% (27%) for high-spin (low-spin) prior. 
For context, Fig. 4 also includes two potential NSBH 
candidates discovered previously; GW190814 (Abbott 
et al. 2020d) is a high-S/N event with well-measured 
masses that has a significantly more massive primary 
and a distinctly more massive secondary than either 
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Figure 4. Component masses of GW200105 (red) 


and GW200115 (blue), represented by their two- and one- 
dimensional posterior distributions. Colored shading and 
solid curves indicate the high-spin prior, whereas dashed 
curves represent the low-spin prior. The contours in the main 
panel, as well as the vertical and horizontal lines in the top 
and right panels, respectively, indicate the 9096 credible in- 
tervals. Also shown in gray are two possible NSBH events, 
GW190814 and the marginal candidate GW190426.152155, 
the latter overlapping GW200115. Lines of constant mass 
ratio are indicated in dashed gray. The green shaded curves 
in the right panel represent the one-dimensional probability 
densities for two estimates of the maximum NS mass, based 
on analyses of nonrotating NSs (Mmax,rov; Landry et al. 
2020, 2021) and Galactic NSs (Mmax,ens; Farr & Chatzi- 
ioannou 2020). 


GW200105 or GW200115, and the marginal candidate 
GW190426.152155 (Abbott et al. 2021c), has (if of as- 
trophysical origin) m1—mz contours that overlap those 
of GW200115. The masses of GW190426.152155 are 
less constrained than those of GW200115 due to its 
smaller S/N. To highlight how the secondary masses of 
GW200105 and GW200115 compare to the maximum 
NS mass, we also show two estimates of the maximum 
NS mass based on an analyses of non-rotating (Landry 
et al. 2020) and Galactic (Farr & Chatziioannou 2020) 
NSs. 

The secondary masses are consistent with the maxi- 
mum NS mass, which we quantify in §5.2. 


4.2. Sky location, distance, and inclination 


We localize GW200105’s source to a sky area of 
7200 deg? (90% credible region). The large sky area 
arises due to the absence of data from LIGO Hanford. 
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Figure 5. Two- and one-dimensional posterior distributions 
for distance Dy, and inclinaton ፀ/ አ. The solid (dashed) lines 
indicate the high-spin (low-spin) prior analysis, and the shad- 
ing indicates the posterior probability of the high-spin prior 
analysis. The contours in the main panel and the horizontal 
lines in the right panel indicate 90% credible intervals. 


The luminosity distance of the source is found to be 
Dr = 280” 116 Mpc. For the second event, GW200115, 
we localize its source to be within 600deg?. It is bet- 
ter localized than GW200105 by an order of magni- 
tude, since GW200115 was observed with three detec- 
tors. We find the luminosity distance of the source to 
be Dy, = 300* 139 Mpc. 

'The luminosity distance is degenerate with the incli- 
nation angle 0; w between the line of sight and the bina- 
ries’ total angular momentum vector (Cutler & Flana- 
gan 1994; Nissanke et al. 2010). Inclination 0yw = 0 
indicates that the angular momentum vector points to- 
ward Earth. The posterior distribution of the inclination 
angle is bimodal and strongly correlated with luminos- 
ity distance, as shown in Figure 5. The inclination mea- 
surement for GW200105 equally favors orbits that are 
either oriented toward or away from the line of sight. 
In contrast, GW200115 shows modest preference for an 
orientation 0; w < 7/2. 


4.3. Spins 


The angular momentum vector S; of each compact 
object is related to its dimensionless spin vector Y; = 
cS;/(Gm2). Its magnitude y; = |YX;| is bounded by 1. 
For GW200105, we infer x1 = 0.08*0 05, which is con- 
sistent with zero. For GW200115, the spin magnitude 


is not as tightly constrained, y1 = (E cur but is also 
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Two-dimensional posterior probability for the spin-tilt angle and spin magnitude for the primary objects 


(left hemispheres) and secondary objects (right hemispheres) for both events. Spin-tilt angles of 0° (180°) correspond to 
spins aligned (antialigned) with the orbital anglular momentum. The color indicates the posterior probability per pixel of 
the high-spin prior analysis. For comparison with the low-spin analysis, the solid (dashed) lines indicate the 90% credible 
regions of the high-spin (low-spin) prior analyses. The tiles are constructed linearly in spin magnitude and the cosine of the 
tilt angles such that each tile contains an identical prior probability. The probabilities are marginalized over the azimuthal angles. 


consistent with zero. The spin of the secondary for both 
events is unconstrained. 

One of the best-constrained spin parameters is the 
effective inspiral spin parameter xeg (Damour 2001; 
Racine 2008; Santamaría et al. 2010; Ajith et al. 2014; 
Vitale et al. 2017a). It encodes information about the 
binaries’ spin components parallel to the orbital angular 
momentum, Xeg — (Yi ፦ 3272) - Ê, where L is the 
unit vector along the orbital angular momentum. 

For GW200105, Xe = —0.01*0:12 and we find the ef 
fective inspiral spin parameter to be strongly peaked 
about zero, with roughly equal support for being ei- 
ther positive or negative. For GW200115, we find mod- 
est support for negative effective inspiral spin: Xef = 
—0.19*022. Negative values of ye indicate binaries 
with at least one spin component negatively aligned 
with respect to the orbital angular momentum, i.e. 
Xiz = Xi: È < 0. We find xi, = —0.19*020, and a 
probability of 88% that x1,; < 0. 

The joint posterior probability of the dimensionless 
spin angular momentum magnitude and tilt angle for 
both components of both events is shown in Fig. 6. The 
tilt angle with respect to the orbital angular momentum 
is defined as arccos (Lgi). Deviations from uniform 
shading indicate a spin orientation measurement. The 
spin orientation of the primary of GW200105 is uncon- 
strained, whereas the orientation of GW200115 shows 
support for negatively aligned primary spin. 


Orbital precession is caused by a spin component in 
the orbital plane of a binary (Apostolatos et al. 1994), 
which we parameterize using the effective precession spin 
parameter 0 < Xp < 1 (Schmidt et al. 2015). We infer 
Xp = 0.0970 bf for GW200105 and xy = 0.217930 for 
GW200115. To assess the significance of a measurement 
of precession, we compute a Bayes factor between a pre- 
cessing and nonprecessing signal model and the preces- 
sion S/N pp (Fairhurst et al. 2020a,b). For GW200105, 
we find a log Bayes factor in favor of spin precession of 
10816 B = —0.24 and precession S/N pp = 0.741 6'37. For 
GW200115, 10816 B = —0.12 and pp = 0.974525. For 
both events and both diagnostics, this indicates incon- 
clusive evidence of precession. This result is expected 
given the S/Ns and inferred inclination angles of the 
binaries (Vitale et al. 2014; Green et al. 2021; Pratten 
et al. 2020b). 

Low values of the primary mass of GW200115 (m4 S 
5 Mo) are strongly correlated with negative values of 
the primary parallel spin component x1,;, as shown in 
Fig. 7. The astrophysical implications of the mass and 
spin correlation are discussed in 86. Figure 7 also shows 
the in-plane spin component x, , which is peaked about 
zero. The lack of conclusive evidence for spin precession 
in GW200115 is consistent with the measurement of x | . 
Apparent differences between the probability density of 
the primary spin in Fig. 6 and the posteriors of x11-X1,z 
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Figure 7. Properties of the primary component of 


GW200115. The corner plot shows the one-dimensional (di- 
agonal) and two-dimensional (off-diagonal) marginal poste- 
rior distributions for the primary’s mass and perpendicular 
and parallel spin components. The shading indicates the 
posterior probability of the high-spin prior analysis. The 
solid (dashed) lines indicate the 50% and 90% credible re- 
gions of the high-spin (low-spin) prior analyses. The vertical 
lines indicate the 90% credible intervals for the analyses with 
high-spin (solid lines) and low-spin (dashed lines) prior. 


in Fig. 7 arise from different choices in visualizing the 
spin orientation posteriors. 


4.4. Remnant properties 


Under the hypothesis of NSBH coalescence for the 
two events, estimates for the final mass and final spin 
of the remnant BH can be made using the models of 
Zappa et al. (2019). We use samples obtained by com- 
bining those from Phenom NSBH and EOB NSBH. 
For GW200105, the remnant mass and spin are Mẹ = 
10.4427 and ye = 0.434903, while for GW200115, Mr = 
7.8114 and yr = 0.384035. We do not investigate any 
post-merger GW signals. The S/Ns of GW200105 and 
GW200115 are around a factor of 3 less than that of 
GW170817, for which there was no evidence of GWs af- 
ter the merger (Abbott et al. 2017b). In the absence of 
tidal disruption, the postmerger signals of GW200105 
and GW200115 would likely resemble a BH ringdown 
(Foucart et al. 2013). The GW signal associated with 
such ringdowns would appear well outside of LIGO’s and 
Virgo’s sensitive bandwidth given the remnant masses 
and spins of the systems (Sarin & Lasky 2021). 


4.5. Tests of general relativity and higher-order GW 
multipole moments 


Results from parameterized tests of general relativ- 
ity (GR; Blanchet & Sathyaprakash 1995; Yunes & Pre- 
torius 2009; Mishra et al. 2010; Li et al. 2012a,b; Agathos 
et al. 2014; Meidam et al. 2018; Abbott et al. 2019b), 
show that GW200105 and GW200115 have too low an 
S/N to allow for tighter constraints than those already 
presented in Abbott et al. (2021d). Within their mea- 
surement uncertainties, our results do not show statis- 
tically significant evidence for deviations from the pre- 
diction of GR. 

To quantify the evidence for higher-order GW multi- 
pole moments, we calculate the orthogonal optimal S/N 
pù, for the subdominant multipole moments (Abbott 
et al. 2020c,d; Mills & Fairhurst 2021). We find the 
(£,|m|) = (3,3) to be the loudest subdominant multi- 
pole moment, as expected for binaries with asymmet- 
ric masses. Using the Phenom HM waveform model, 
we infer pł = 1.70*092 (1.70*995) for GW200105 and 
33 = 0.91* 005 (0.86*0:05) for GW200115 with the low 
(high) spin prior. In Gaussian noise, the median of p33 is 
approximately chi-distributed with two degrees of free- 
dom, and values greater than 2.1 indicate significant 
higher-order multipole content. The measured p33 are 
therefore consistent with Gaussian noise, as expected for 
the majority of NSBHs at these S/Ns, except for those 
viewed close to edge-on. 


4.6. Waveform systematics 


Our primary results are obtained using precessing 
BBH models with higher-order multipole moments, Phe- 
nom PHM and EOBNR PHM. We now justify this 
choice by investigating potential systematic uncertain- 
ties due to our waveform choice. 

First, we investigate the agreement between inde- 
pendent waveform models that incorporate identical 
physics. Figure 8 shows the two-dimensional m2—Veg 
posteriors for both events obtained using a variety of 
NSBH and aligned-spin BBH models. Because some 
NSBH models only cover xı < 0.5, we restrict the prior 
range of all models to xı < 0.5 for consistency. 

'The main panels of Fig. 8 are dominated by a corre- 
lation of the effective inspiral spin parameter xeg with 
the secondary mass mz (Cutler & Flanagan 1994; Ng 
et al. 2018a). Both NSBH models (Phenom NSBH and 
EOBNR NSBH) give consistent results with each other, 
as do both BBH models (Phenom and EOBNR), both 
with and without higher-order multipole moments, with 
the most notable difference being that EOBNR HM 
yields tighter posteriors than Phenom HM. This demon- 
strates that waveform models including the same physics 
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Figure 8. Comparison of two-dimensional m2-Xef posteriors for the two events reported here, using various NSBH and BBH 
signal models. The vertical dashed lines indicate several mass-ratio references mapped to m» for the median estimate of the 


chirp masses of GW200105 and GW200115. 


give comparable results, but more studies are warranted 
to improve the understanding of the BBH waveform 
models in the NSBH region of parameter space. While 
not shown in Fig. 8, we also find good agreement be- 
tween the primary precessing BBH waveform models; 
see §4. 

Second, comparing the NSBH models with the BBH 
models without higher-order multipole moments (Phe- 
nom and EOBNR)), the NSBH models recover similar 
posterior contours in the mo-xegr plane. This is ex- 
pected given the asymmetric mass ratio and low S/N of 
these NSBH observations; see, e.g. Huang et al. (2021) 
for a demonstration that higher S/Ns would be needed 
to see notable systematic effects. We observe differences 
at the extreme ends of the m»-xeg contours (1.6. at the 
smallest and largest values of m2). The construction 
of the NSBH waveform models used here did not rely 
on numerical relativity results at mass ratios q S 1/8, 
nor did they include simulations with x1: < 0 or NS 
masses m3 > 1.4.7, (Matas et al. 2020; Thompson 
et al. 2020). Therefore, some differences should be ex- 
pected, especially for large mz in GW200105. Further- 
more, for GW200105, the tails of the mo-xeg distribu- 
tion for Phenom NSBH and EOB NSBH at high mə are 
also impacted by the inability of the data to constrain 
the tidal deformability. Hence, the posterior samples 
include combinations of high mz with large Ag, despite 
such combinations being unphysical. This effect is not 
apparent for GW200115 because of its smaller secondary 


mass. The isolated islands of probability in the extreme 
tails of the distributions are due to sampling noise. 
Last, when adding the extra physical content 
of higher-order multipole moments in BBH models 
(through Phenom HM and EOBNR HM), the extreme 
ends of the ma-xeg contours are excluded, while the 
bulk of the distributions are consistent with the posteri- 
ors obtained with the NSBH models. In summary, these 
comparisons indicate that (i) waveform models includ- 
ing the same physics give comparable results; (ii) going 
from NSBH models to comparable BBH models changes 
the results only marginally, i.e. any effects of tides are 
small; and (iii) inclusion of higher-order multipole mo- 
ments changes the posterior contours more substantially 
than inclusion of tides. We conclude that the the inclu- 
sion of precession and higher-order multipole moments 
afforded by the BBH waveform models is more impor- 
tant than the impact of tides in the NSBH models. 


5. NATURE OF THE SECONDARY COMPONENTS 


In this section, we describe the investigations to es- 
tablish the nature of the secondary objects. In 85.1, we 
look for imprints of tidal deformations of the secondaries 
and conclude that the masses, spins, distances and S/Ns 
of the detections make definitive identifications of NSs 
unlikely, both in GW and EM measurements. However, 
in 85.2, we show that the posterior distributions of the 
secondary masses agree with those of known NSs. 
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5.1. Tidal deformability and tidal disruption 


The tidal deformability of NSs is imprinted in the GW 
signal, and is investigated using the parameter estima- 
tion techniques of $4. In contrast, BHs have zero tidal 
deformability (Binnington & Poisson 2009; Damour & 
Nagar 2009; Chia 2020; Charalambous et al. 2021; Le 
Tiec & Casals 2021). We infer the tidal deformabil- 
ity Ag of the NSs in GW200105 and GW200115 using 
the NSBH waveform models that include tides. We find 
that the tidal deformabilities are uninformative relative 
to a uniform prior in A» € [0, 5000]. This measurement 
cannot establish the presence of NSs, which is expected 
given the mass ratios and the S/N of the detections (Fou- 
cart et al. 2013; Kumar et al. 2017; Fasano et al. 2020; 
Huang et al. 2020). 

Toward the end of the inspiral, the BH may tidally 
disrupt the NS and form an accretion disk (Pannar- 
ale 2013; Foucart et al. 2018). This is hypothesized to 
drive a relativistic jet (Pannarale et al. 2011; Pascha- 
lidis et al. 2015). Given the mass ratios for both events 
and the aligned spins x1,: of their primaries (near zero 
for GW200105, probably negative for GW200115), we 
do not expect tidal disruption to occur, which would re- 
quire more equal masses or more positive x1, (Rantsiou 
et al. 2008; Shibata & Taniguchi 2008; Etienne et al. 
2009; Foucart et al. 2011; Kyutoku et al. 2011; Pannar- 
ale 2013; Foucart et al. 2018). 

To quantitatively confirm this expectation, we use 
the spectral representation of equations of state from 
Abbott et al. (2018a), which uses an SLY low-density 
crust model (Douchin & Haensel 2001) and parame- 
terizes the adiabatic index into a polynomial in the 
logarithm of the pressure (Lindblom 2010; Lindblom 
& Indik 2012, 2014). Following Stachie et al. (2021), 
we marginalize the parameter estimation samples from 
the NSBH analyses over these equations of state. For 
a fixed equation of state, we compute the maximum 
Tolman-Oppenheimer-Volkov (TOV) mass, allowing us 
to infer the nature (NS or BH) of the lighter binary 
compact object, as well as its radius Ro, compactness 
Cz = Gmzs/(Ra5c?) and baryon mass. Based on these, 
we define a total ejecta mass mej (Fernández et al. 2019) 
as the sum of dynamical ejecta (Krüger & Foucart 2020) 
and 1596 of the mass of disk winds (Foucart et al. 2018). 
For both events, we find that mej < 10-9 Mc for 99% of 
the samples. 

'The absence of ejecta is compatible with the lack of 
observed EM counterparts. However, given the large 
distances of the mergers (~ 300 Mpc) and the large un- 
certainties of their sky localization, EM emission would 
have been difficult to detect and associate with these 
GW events. 


Table 3. Probability that the secondary mass is below the 
maximum NS mass Mmax for each event, given different spin 
assumptions and different choices for the maximum NS mass. 
The values shown use a flat prior in ma; alternative, astro- 
physically motivated mass priors, can cause the estimates to 
vary by up to 1196 across our chosen models. 


p(ma « Mmax) 
spin prior choice of Mmax GW200105 GW200115 


22] < 0.05 Mmax,TOV 96% 98% 
|x2| < 0.99 Mmax(X2) 94% 95% 
Ix2] < 0.99 — Mmax,aNs 93% 96% 


Estimating the impact of nonlinear p-g tidal coupling 
(Abbott et al. 2019c), we find that it would produce 
a relative frequency-domain phase shift for GW200105 
(GW200115) of approximately 134 (38) times smaller 
than the equivalent phase shift for GW170817. This 
strongly reduced effect is caused by the larger chirp 
masses, more asymmetric mass ratios, and the presence 
of only a single NS. Since p-g effects were not detected 
for GW170817 (Abbott et al. 2019c), they will be unob- 
servable within the new NSBH systems. 


5.2. Consistency of component masses with the NS 
maaimum mass 


Even without definite identification of matter signa- 
tures in the signals, we can compare the observed m» 
for GW200105 and GW200115 with the maximum NS 
mass, Mmax. The existence of massive pulsars (Anto- 
niadis et al. 2013; Cromartie et al. 2019; Fonseca et al. 
2021) places a lower bound of Mmax Z 2Mo on the 
maximum NS mass. Studies of GW170817's remnant 
typically suggest that the maximum mass of a nonrotat- 
ing NS—the TOV mass—is Mmax,rov S 2.3 Mo (e.g., 
Shibata et al. 2019; Abbott et al. 2020b; Nathanail et al. 
2021). However, rapid rotation could support a larger 
Mmax. Given the considerable uncertainty in Mmax, we 
examine three different scenarios. Following Essick & 
Landry (2020), we compute for each scenario the prob- 
ability p(mo < Mmax) that the secondary mass is below 
the maximum NS mass by marginalizing over the uncer- 
tainty in Mmax and the uncertainty of our mo measure- 
ment. 

Supposing that the secondaries are slowly spinning, 
we consider in the first row of Table 3 an estimate of 
Mmax,rov from a nonparametric astrophysical inference 
of the equation of state (Landry et al. 2020), which pre- 
dicts Mmax,rov = 2.22*030 Mo and is shown in Fig. 4. 
We then relax the low-spin assumption, estimating in 
the second row the maximum rotationally supported 


mass Mmax(X2), and the breakup spin Xmax with the 
universal relations from Breu & Rezzolla (2016). In this 
scenario, we require m3 < Minax(x2) and %2 < Xmax 
for consistency with an NS. Finally, in the third row, 
we consider a parametric fit to the entire distribution 
of observed Galactic NSs, including rapidly rotating 
pulsars (Farr & Chatziioannou 2020), which predicts 
Mmax ans = 2.25*0 11 Mo, and is shown in Fig. 4. This 
scenario accounts for the possibility that the maximum 
mass in the NS population is limited by the astrophysi- 
cal processes that form compact binaries. Assuming low 
spin (first row), we find probabilities of 96% and 98% 
that the secondaries in GW200105 and GW200115, re- 
spectively, are consistent with an NS assuming a uniform 
prior in m2. The possibility of large secondary spin re- 
duces these probabilities by up to 396 (second and third 
rows). 

So far, this analysis has assumed priors that are uni- 
form in component masses. However, there is consid- 
erable uncertainty in the astrophysical mass priors of 
such systems, and different prior assumptions can af- 
fect the component mass posteriors for detections with 
moderate S/N. To illustrate the impact of population 
assumptions, we consider three alternative priors: one 
based on Salpeter mass distributions, p(m) ~ m ?? 
(Salpeter 1955), independently for each component; one 
based on an extrapolation of the BBH mass model BRo- 
KEN PowER LAW from Abbott et al. (2021b) down to 
0.5 Mo for both components; and another based on a 
similar extrapolation of the POWER LAW + PEAK BBH 
mass model from the same study. We marginalize over 
the uncertainties in the latter two models, which are fit 
to the BBH population from Abbott et al. (2021c), in- 
cluding the outlier event GW190814 with a secondary 
component mass below 3 Mg (Abbott et al. 2020d). 

These different mass priors change the numbers in 
Table 3 by at most 1196, with the smallest values for 
GW200105 and GW200115 being 89% and 87%, respec- 
tively. The decrease is due to the three priors assigning 
more probability density to equal-mass systems, thus fa- 
voring larger ma. Thus, the secondaries of both systems 
are consistent with NSs based on our assumptions about 
the equation of state and the mass distribution. 

However, consistency with the maximum NS mass 
does not exclude the possibility that the secondaries 
could be BHs or exotic compact objects, if such objects 
also exist within the NS mass range. For instance, mod- 
els of primordial BHs predict a peak in the primordial 
BH mass function at ~ 1 Mọ (Carr et al. 2021). These 
models also predict that primordial BHs may form co- 
alescing binaries at mass ratios comparable to those re- 
ported here. 
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Figure 9. Inferred probability densities for the NSBH 


merger rate. Green line: rate assuming one count each from 
an GW200105 and GW200115-like NSBH population. Black 
line: rate for a broad NSBH population with a low threshold 
that accounts for marginal triggers. The short vertical lines 
indicate the 90% credible intervals. 


6. ASTROPHYSICAL IMPLICATIONS 


The first confident observations of NSBH binaries en- 
able us to study this novel type of astrophysical system 
in entirely new ways. We pursue three different avenues 
in this section. First, we infer the merger rate of NSBH 
binaries in the local universe. We then place the inferred 
source properties and merger rate in the context of mod- 
els of NSBH formation channels and previous EM and 
GW observations of BHs and NSs. Finally, we investi- 
gate to what extent the events reported here can serve 
to measure the Hubble constant and whether lensing of 
GWs may have played a role in the observations. 


6.1. Merger rate density 


We infer the NSBH merger rate density with our ob- 
servations using two different approaches. 

In the first approach, we consider only GW200105 
and GW200115. Following the method of Kim et al. 
(2003) as previously used in, e.g., Abbott et al. (20166, 
2020d), we calculate an event-based merger rate assum- 
ing one Poisson-distributed count each from GW200105- 
and GW200115-like populations. We calculate semian- 
alytically the search sensitivities across 01, 02, and 
the first nine months of 03 to NSBH populations cor- 
responding to the mass and spin posteriors for the 
two events. We then calibrate these sensitivities to 
the results of GSTLAL (Cannon et al. 2012; Privitera 
et al. 2014; Messick et al. 2017; Sachdev et al. 2019; 
Hanna et al. 2020) using a broad NSBH-like popula- 
tion and an FAR threshold of 1/(2.8 yr). This FAR 
threshold is chosen to include only GW200105 and 
GW200115 while excluding low-significance triggers like 
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GW190426.152155, though a more conservative thresh- 
old of 1/(100 yr) as used in, e.g., Abbott et al. (20166, 
2020d), changes the estimated sensitivities by less than 
1596. AREE the Poisson Jeffreys prior proportional 
to RC = we find 75300105 = 1625 Gpc^? yr^!, and 
R200115 = 36135 Gpc73 yr-!. The PYCBC detection of 
GW200115, using the same method but with an inde- 
pendent set of injections for calibration (Tiwari 2018), 
yields a consistent 75300115 = 40*22 Gpc^? yr-!. Com- 
bining the likelihoods over R200105 and R200115 accord- 
ing to Kim et al. (2003) and applying the Jeffreys prior 
to the total rate, we find the total event-based NSBH 
merger rate density Rnspy = 45275 Gpc^? yr7!, plot- 
ted in green in Fig. 9. 

The second approach to calculating a merger rate 
takes into account not only GW200105 and GW200115, 
but also less significant search triggers with masses con- 
sistent with the typical range associated with NSBH bi- 
naries. Specifically, we consider triggers across O1, O2 
and the first nine months of O3 with associated com- 
ponent masses m4 € [2.5, 40] Mo and m» € [1,3] Mo, 
ie. broader ranges than the component mass poste- 
riors for GW200105 and GW200115 obtained in 84. 
The cut-off of mz < 3Mg is chosen as a robust up- 
per limit on the maximum NS mass (Rhoades & Ruffini 
1974; Kalogera & Baym 1996). The population is de- 
fined to be uniformly distributed in comoving volume, 
uniform in log component masses in the given ranges, 
with aligned spins with spin magnitudes distributed uni- 
formly in [0,0.95]. We use as input GSTLAL search 
triggers above an S/N threshold such that the number 
of noise triggers exceeds the number of astrophysical sig- 
nals by a factor ~100. Following Farr et al. (2015); Ka- 
padia et al. (2020), we find the resulting joint likelihood 
on Poisson parameters for signals of each astrophysical 
category (Apns, AusBH, Agga), and for terrestrial noise 
triggers Apackground- Here, the BBH and BNS categories 
are defined to include triggers where both component 
masses fall within 5-100 Mo or 1-2.5 Mo, respectively. 
We apply the Jeffreys prior and recover a merger rate 
density Rnspu = AnsBH/(VT)nsBH, where (VT)nsBH 
is the population-averaged sensitive time-volume esti- 
mated using the GSTLAL pipeline and the NSBH pop- 
ulation defined above. This yields an estimated NSBH 
merger rate density of Rnspu = 130166 Gpc-^?yr-!, 
the black line in Fig. 9. 

While this rate is higher than our event-based rate, 
it considers a wider population that includes addi- 
tional triggers; for example, the component masses of 
GW190814 (although the nature of its secondary is un- 
certain) fall within this population. The calculation is 
also based on the component mass values of GSTLAL 


search triggers, adjusted with an analytical model (Fong 
2018) of the response of the template bank to signals 
with a given distribution of true masses. This procedure 
is expected to be less precise than Bayesian parameter 
estimation, which is impractical for the large numbers 
of triggers involved. 

'The merger rate density measured here is consistent 
with the upper bound of 610 Gpc^? yr-! derived from 
the absence of NSBH detections in O1 and O2 (Abbott 
et al. 2019a). Revised merger rate estimates for all CBC 
sources will be provided in a future analysis of the full 
O3 data set. 


6.2. System origins 


To understand the origin of GW200105 and 
GW200115, we compare their observed masses and spins 
with theoretical predictions. Population synthesis stud- 
ies modeling the various formation channels of merging 
compact object binaries distinguish between NSs and 
BHs with a simple mass cut, typically between 2 and 
3 Me. This is consistent with the secondary masses of 
both events being classified as NSs, so for the purposes of 
discussing formation channels for these events and pre- 
dicted rates, we will discuss GW200105 and GW200115 
in the context of them being NSBHs. 


6.2.1. Formation channels 


Formation channels of NSBH can be broadly catego- 
rized as either isolated binary evolution or one of several 
dynamical formation channels (e.g., globular clusters or 
nuclear star clusters). Since isolated binaries form in 
young star clusters and can be influenced by dynamical 
interactions before the cluster dissolves and the binary 
effectively becomes isolated, rates from young star clus- 
ters naturally encompass rates from isolated binaries. 
Predicted rates of NSBH mergers in the local universe 
vary by orders of magnitude across the various forma- 
tion channels. 

Models of the canonical isolated binary evolution 
channel—in which stellar progenitors evolve together, 
shedding orbital angular momentum through phases of 
stable and/or unstable mass transfer prior to compact 
object formation—predict NSBH merger rate densities 
around 0.1-800 Gpc^? yr^! (Belezynski et al. 2002; Sip- 
ior & Sigurdsson 2002; Belczynski et al. 2006, 2016; Do- 
minik et al. 2015; Eldridge et al. 2017; Kruckow et al. 
2018; Mapelli & Giacobbo 2018; Neijssel et al. 2019; 
Drozda et al. 2020; Zevin et al. 2020; Broekgaarden et al. 
2021). The high uncertainty is driven by the lack of 
observed NSBHs and the wide range of model assump- 
tions. Merger rates are sensitive to the treatment of 
common envelopes, which may be a necessary evolu- 
tionary phase for producing compact binaries in tight 


orbits capable of merging in a Hubble time (Ivanova 
et al. 2013). They are also sensitive to prescriptions for 
supernova kick magnitudes. While moderate kicks can 
produce eccentric orbits that merge on short timescales, 
high supernova kicks may disrupt the progenitor bina- 
ries and suppress the merger rate (Belczynski et al. 2002; 
Giacobbo & Mapelli 2020; Tang et al. 2020). 

Models of star formation in the dynamical environ- 
ments of young star clusters predict NSBH merger rate 
densities of 0.1-100 Gpc^? yr^! (Fragione & Banerjee 
2020; Hoang et al. 2020; Rastello et al. 2020; Santoliq- 
uido et al. 2020). In this scenario, most systems that 
form merging NSBH (~80%) are ejected without under- 
going dynamical exchanges, proceeding to merge in the 
field. 

Models of dynamical formation channels in denser en- 
vironments typically predict much lower merger rates. 
For instance, in globular clusters and nuclear star clus- 
ters, BHs segregate and dominate the core, where the 
bulk of dynamical interactions occur (Portegies Zwart & 
McMillan 2000; Morscher et al. 2015), so that encounters 
between NSs and BHs are relatively rare, with NSBH 
merger rate densities on the order of 107? Gpc^? yr7! 
(Clausen et al. 2013; Arca Sedda 2020; Ye et al. 2020). 
In disks of active galactic nuclei, the presence of gas 
could possibly increase the NSBH merger rate density 
up to 800 Gpc ^? yr! (McKernan et al. 2020). 

NSBHs may also merge via hierarchical triple inter- 
actions, where inner NSBH binaries are driven to high 
eccentricity by massive tertiary companions and merge 
on rapid timescales (Antonini et al. 2017; Silsbee & 
Tremaine 2017). However, the predicted merger rates 
are negligible unless supernova kicks are assumed to be 
zero (Fragione & Loeb 2019). 

It is likely that a combination of the above channels 
contribute to the astrophysical NSBH merger rate. How- 
ever, the isolated binary evolution, young star cluster, 
and active galactic nuclei channels are capable of indi- 
vidually accounting for the NSBH merger rate estimated 
here. 


6.2.2. Masses 


While there are no observed NSBHs in the Milky 
Way, we can place the component masses of GW200105 
and GW200115 in the context of the observed pop- 
ulation of BH and NS masses, as well as the pre- 
dicted populations of NSBHs. Observations suggest 
that the mass distribution of the Galactic population 
of NSs peaks around 1.33 Mo, with a secondary peak 
around 1.9 Mc (Antoniadis et al. 2016; Alsing et al. 
2018). The secondary mass observed in GW200115 and 
marginal event GW190426 152155 are consistent with 
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the population peaking at 1.33 Mo, while the secondary 
observed in GW200105 (c 1.9 Mo) and the primary 
component from BNS merger GW190425 (m; = 1.60- 
1.87 Me; Abbott et al. 20208) are consistent with the 
high-mass population. However, a rigorous association 
of the events with different components of the NS pop- 
ulation would require a thorough population analysis. 
Radio observations of BNS systems do not find such 
massive NSs, leading to speculation as to the origin of 
GW190425 (Romero-Shaw et al. 2020; Safarzadeh et al. 
2020; Galaudage et al. 2021; Mandel et al. 2021). Stel- 
lar metallicities in the Milky Way are not representative 
of all populations of GW sources (O'Shaughnessy et al. 
2008, 2010; Belczynski et al. 2010; Eldridge et al. 2017; 
Neijssel et al. 2019). 

The BH masses observed in GW200105 and 
GW200115 (8.9*12 Mo and 5.7*11 Mo, respectively) 
are in line with predictions from population synthesis 
models for NSBH mergers from isolated binary evolu- 
tion and young star clusters. In NSBHs, current binary 
evolution models do not predict BH masses above ~ 10 
Mo (Eldridge et al. 2017; Kruckow et al. 2018; Mapelli & 
Giacobbo 2018; Broekgaarden et al. 2021), while Popu- 
lation III NSBHs (Kinugawa et al. 2017) and dynamical 
interactions in low-metallicity young star clusters allow 
for higher BH masses (Ziosi et al. 2014; Rastello et al. 
2020). 

Electromagnetic observations of X-ray binaries have 
not uncovered BHs between 3 and 5 Mo, leading to spec- 
ulation about a mass gap (Ozel et al. 2010; Farr et al. 
2011; Kreidberg et al. 2012; Miller & Miller 2014). Anal- 
ysis of GWTC-2 has also found evidence for a gap or 
dip in the BH mass spectrum between ~2.6 and 4 Mo 
(Fishbach et al. 2020). For GW200115, we find nonneg- 
ligible support for the primary lying in this mass gap, 
with p(3 Mo < mi < 5 Mc) = 30% (27%) under the 
high-spin (low-spin) priors. This low-mass region is cor- 
related with negative values of the parallel component 
of the primary spin; see $6.2.3 below. 

In summary, the masses inferred for GW200105 and 
GW200115 are consistent with expectations for NSBHs; 
their primary masses are in agreement with predictions 
for BH masses in population synthesis models of the 
dominant formation scenarios. Meanwhile, their sec- 
ondary masses are compatible with the observed pop- 
ulation of Galactic NSs, as well as the masses inferred 
from GW observations of BNS mergers. 


6.2.3. Spins 


Spin information encoded in GWs from binaries are 
a probe of the their evolutionary history (Farr et al. 
2017, 2018; Stevenson et al. 2017; Talbot & 'Thrane 
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2017; Vitale et al. 2017b; Wysocki et al. 2019). The 
BHs in binaries are expected to exhibit a range of spin 
magnitudes and orientations, depending on how they 
formed (Kalogera 2000; Mandel & O'Shaughnessy 2010; 
Rodriguez et al. 2016; Liu & Lai 2017, 2018; Antonini 
et al. 2018; Kruckow et al. 2018; Bavera et al. 2020; 
Chattopadhyay et al. 2020). The highest dimensionless 
NS spin implied by pulsar-timing observations of bina- 
ries that merge within a Hubble time is 0.04 (Stovall 
et al. 2018; Zhu et al. 2018). 

While the secondary spins of both events reported 
here are poorly constrained due to the unequal masses, 
the primary spins of GW200105 and GW200115 can 
be placed in the context of predictions for BH spins 
from models of stellar and dynamical evolution and 
EM observations of NSBH progenitors. As can be seen 
in Fig. 6 and Fig. 7, the primary spins of GW200105 
and GW200115 are consistent with zero (0.08*022 and 
033025. respectively), but moderate values of spin are 
not ruled out. The primary in GW200115 may even have 
relatively high spin, with a 9096 upper limit of 0.72. Sev- 
eral studies of the observed population of high-mass X- 
ray binaries (Liu et al. 2008; Gou et al. 2009, 2014; Zhao 
et al. 2021) find that the BHs have large spins (Valsec- 
chi et al. 2010; Wong et al. 2012; Qin et al. 2019; Zhao 
et al. 2021). Given the short lifetime of the secondary, 
mass transfer is argued to be insufficient to generate BHs 
with such high spins, implying that the BHs were born 
with high spins. Belczynski et al. (2011) found that one 
such high-mass X-ray binary, Cygnus X-1, is expected 
to form an NSBH with a BH that carries near-maximal 
spin, although it would not merge within a Hubble time. 
However, following revised estimates of the component 
masses of Cygnus X-1 (Miller-Jones et al. 2021), Nei- 
jssel et al. (2021) found that it will most likely form a 
BBH. Meanwhile, analyses of GWTC-1 and GWTC-2 
have found evidence for BH spin (Abbott et al. 2016f; 
Vitale et al. 2017; Chatziioannou et al. 2019; Kimball 
et al. 2020; Abbott et al. 2021b), though they do not 
determine whether those BHs may have been formed 
with that spin. Altogether, these EM and GW observa- 
tions of compact binaries and their progenitors suggest 
a range of BH natal spins in NSBH binaries. 

Along with their magnitudes, the alignments of com- 
ponent spins with the overall binary orbital angular mo- 
mentum are of astrophysical interest. In particular, we 
find evidence that the primary BH spin in GW200115 
is negatively aligned with respect to the orbital angular 
momentum axis, with p(xi,; < 0) = 88% (87%) under 
the high-spin (low-spin) prior and with the more neg- 
ative values of x1,, correlated with smaller mi. This 
negative alignment is consistent with dynamical forma- 


tion channels, which typically form binaries with ran- 
dom spin orientations (Rodriguez et al. 2016), but the 
predicted rates from these channels, discussed in 86.2, 
are small. Binaries born in isolation are expected to 
form with only small misalignments (ፎ 30°; Kalogera 
2000), though they may become misaligned by super- 
nova kicks at compact object formation (Rodriguez et al. 
2016; Gerosa et al. 2018; Wysocki et al. 2018), and 
possibly during subsequent evolution via mass transfer 
(Stegmann & Antonini 2021). Meanwhile, NSBH pro- 
genitor binaries originating in young star clusters can 
be perturbed via close dynamical encounters before be- 
ing ejected into the field. Therefore, a misaligned spin 
in the primary of GW200115 would not necessarily rule 
out any of the plausible NSBH formation channels. 


6.3. Cosmology and lensing 


Gravitational wave sources are standard sirens, pro- 
viding a direct measurement of their luminosity dis- 
tance (Schutz 1986; Holz & Hughes 2005), and they can 
be used to measure the Hubble constant (Abbott et al. 
2017c; Fishbach et al. 2019; Abbott et al. 2021a). Due to 
the lack of a confirmed EM counterpart and large num- 
bers of galaxies inside the localization volumes of each of 
the two events, we do not obtain any informative bounds 
on Ho from these observations. 

The detections of GW200105 and GW200115 are sep- 
arated by only ^10 days. One explanation for the small 
time-delay could be that the two events are created by 
gravitational lensing by a galaxy (Haris et al. 2018). 
Gravitational lensing is unlikely even a priori (Smith 
et al. 2017; Ng et al. 2018b), and the nonoverlapping 
mass posteriors (Fig. 4) further exclude it as a possible 
explanation (Haris et al. 2018). While GW200115 and 
GW190426_152155 exhibit agreement in their source 
mass posteriors, their sky localization areas do not over- 
lap, and their detector-frame (redshifted) chirp masses 
show only marginal overlap (Abbott et al. 2021c), ruling 
out lensing as a possible explanation. 


7. CONCLUSIONS 


During its third observing run, the LIGO-Virgo GW 
detector network observed GW200105 and GW200115, 
two GW events consistent with NSBH coalescences. 
GW200105 is effectively a single-detector event observed 
in LIGO Livingston with an S/N of 13.9. It clearly 
stands apart from all recorded noise transients, but its 
statistical confidence is difficult to establish. GW200115 
was observed in coincidence by the network with an S/N 
of 11.6 and FAR of <1/(1 x 10° yr). 

The source component masses of GW200105 and 
GW200115 make it likely that these events originated 


from NSBH coalescences. Their primary masses are 
found to be m, = 8.9*12 Mo and mı = 5.7154 Mo, 
which are consistent with predictions of BH masses in 
population synthesis models for NSBHs. Their sec- 
ondary masses, inferred to be mg = 1.9733 Me and 
m = Lan: Mo, respectively, are consistent with the 
observed NS mass distribution in the Milky Way, as well 
as population synthesis predictions for secondary masses 
in merging NSBHs. 

We find no evidence of measurable tides or tidal dis- 
ruption for either of the two signals, and no EM coun- 
terparts to either detection have been identified. As 
such, there is no direct evidence that the secondaries 
are NSs, and our observations are consistent with either 
event being a BBH merger. However, the absence of 
tidal measurements and EM counterparts is to be ex- 
pected given the properties and distances of the two 
events. Moreover, the comparisons of the secondary 
masses to the maximum allowed NS mass yield a proba- 
bility p(m» € Mmax) of 8996-0696 and 8796-9896 for the 
secondaries in GW200105 and GW200115, respectively, 
being compatible with NSs (see 85.2). 

The effective inspiral spin parameter of GW200105 is 
strongly peaked around zero: xeg — —0.017012. For 
the second event, GW200115, the effective inspiral spin 
parameter is inferred to be Xef = UNDIS For 
GW200115, the spin component parallel to the orbital 
angular momentum of the primary is x1,; = -0197$ 27; 
and we find support for negatively aligned primary spin 
(Xı,z < 0) at 88% probability. More negative values of 
X1,z in GW200115 are correlated with particularly small 
primary masses reaching into the lower mass gap. We 
find p(3 Mo < mı «5 Mo) = 3096 (27%) under the high- 
spin (low-spin) parameter estimation priors. We find no 
conclusive evidence for spin-induced orbital precession 
in either system. 

We estimate the merger rate density of NSBH bina- 
ries with two approaches. Assuming that GW200105 
and GW200115 are representative of the entire NSBH 
population, we find Ruspu = 45* 25 Gpc^? yr-!. Con- 
versely, assuming a broader range of allowed primary 
and secondary masses, and considering all triggers in 
03, we find Rwsagg = 130766. Gpc73yr7!. These 
are the first direct measurements of the NSBH merger 
rate. Both estimates are broadly consistent with the 
rate predicted from NSBH formation in isolated bina- 
ries or young star clusters. Formation channels in dense 
star clusters (globular or nuclear) and in triples predict 
lower rates than those inferred from the two events and 
are unlikely to be the dominant NSBH formation chan- 
nels. 
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The observations of GW200105 and GW200115 are 
consistent with predictions for merging NSBHs and ob- 
servations of BHs and NSs in the Milky Way. Given their 
significantly unequal component masses, future observa- 
tions of NSBH systems will provide new opportunities to 
study matter under extreme conditions, including tidal 
disruption, and search for potential deviations from GR. 
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sequent significance evaluation were performed with 
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2017; Sachdev et al. 2019; Hanna et al. 2020), built 
on the LALSUITE software library (LIGO Scientific 
Collaboration 2018), and with the PyCBC (Usman 
et al. 2016; Nitz et al. 2018, 2019) and MBTAOn- 
LINE (Adams et al. 2016) packages. Parameter estima- 
tion was performed with the LALINFERENCE (Veitch 
et al. 2015) and LALSIMULATION libraries within LAL- 
SUITE (LIGO Scientific Collaboration 2018), as well as 
the BILBY and PBILBY Libraries (Ashton et al. 2019; 
Smith et al. 2020) and the DYNESTY nested sampling 
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tra were obtained using BAYESWAVE (Cornish & Lit- 
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ter estimation results were visualized and shared with 
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2020). Plots were prepared with Matplotlib (Hunter 
2007). Figure 1 was generated using GWpy (Macleod 
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2013; Price-Whelan et al. 2018) and ligo.skymap 
(https: //Iscsoft.docs.ligo.org/ligo.skymap). 

We would like to thank all of the essential workers who 
put their health at risk during the COVID-19 pandemic, 
without whom we would not have been able to complete 
this work. 
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